clear all
close all
hold off
clc
bval = 1 ;


clearvars -except bval
eval(['load(''data/dataformatlab_new_2b' int2str(bval) '.mat'')']) 
run set_parameters.m 

indvar = linspace(1,N*2,N*2)-1 ;
[xhat4,fval4,exitflag4]         = fminsearch(@(x) adcurve5(x,indvar(1:N),CumulativeMeanAcceptedCompMen(1:N)./CumulativeMeanAcceptedCompMen(1)),[.7 .1],options) ;
cumcomp_m_dat                   = CumulativeMeanAcceptedCompMen(1).*(xhat4(1)+(1-xhat4(1)).*exp(-xhat4(2).*(indvar))) ;
plot(indvar(1:N),cumcomp_m_dat(1:N),'r',indvar(1:N),CumulativeMeanAcceptedCompMen(1:N),'bo')

indvar5                         = linspace(1,N*2,N*2)-2 ;
[xhat6]                         = fminsearch(@(x) adcurve5(x,indvar5(5:N),CumulativeMeanAcceptedCompWomen(5:N)./CumulativeMeanAcceptedCompWomen(5)),[.7 .1],options) ;
cumcomp_f_dat                   = CumulativeMeanAcceptedCompWomen(5).*(xhat6(1)+(1-xhat6(1)).*exp(-xhat6(2).*(indvar5))) ;
plot(indvar5(1:N),cumcomp_f_dat(1:N),'r',indvar5(1:N),CumulativeMeanAcceptedCompWomen(1:N),'bo',...
    indvar(1:N),cumcomp_m_dat(1:N),'r',indvar(1:N),CumulativeMeanAcceptedCompMen(1:N),'bo')

gendergap_dat   = cumcomp_m_dat./cumcomp_f_dat ;
plot(gendergap_dat)
[xhat8,fval8,exitflag8]         = fminsearch(@(x) adcurve6(x,indvar(1:N),CumulativeProportionMen(1:N)),[0 CumulativeProportionMen(N) .01 CumulativeProportionMen(1) 10 1 1],options) ;
[cumsum_m_dat]                  = min(max((xhat8(1)+(xhat8(2)-xhat8(1))./(xhat8(3)+xhat8(4).*exp(-xhat8(5).*(indvar'-xhat8(6)))).^(1./xhat8(7))),0),1) ;
plot(indvar(1:N),cumsum_m_dat(1:N),'r',indvar(1:N),CumulativeProportionMen(1:N),'b')

[xhat10,fval10,exitflag10]      = fminsearch(@(x) adcurve6(x,indvar(1:N),CumulativeProportionWomen(1:N)),[xhat8],options) ;
[cumsum_f_dat]                  = min(max((xhat10(1)+(xhat10(2)-xhat10(1))./(xhat10(3)+xhat10(4).*exp(-xhat10(5).*(indvar'-xhat10(6)))).^(1./xhat10(7))),0),1) ;
plot(indvar(1:N),cumsum_f_dat(1:N),'r',indvar(1:N),CumulativeProportionWomen(1:N),'b')

[xhat8,fval8,exitflag8]         = fminsearch(@(x) adcurve6(x,indvar(1:20),(AdjustedProbabilityofReceivinganOfferMen(1:20))),[.1 .5  .5 .5 .5 .5 .6],options) ;
[proboffer_m_dat]               = min(max((xhat8(1)+(xhat8(2)-xhat8(1))./(xhat8(3)+xhat8(4).*exp(-xhat8(5).*(indvar'-xhat8(6)))).^(1./xhat8(7))),0),1) ;
plot(indvar(1:N),proboffer_m_dat(1:N),'r',indvar(1:N),(AdjustedProbabilityofReceivinganOfferMen(1:N)),'b')

[xhat8,fval8,exitflag8]         = fminsearch(@(x) adcurve6(x,indvar(1:20),(AdjustedProbabilityofReceivinganOfferWomen(1:20))),[.01 .5  .5 .5 .5 .5 .6],options) ;
[proboffer_f_dat]               = min(max((xhat8(1)+(xhat8(2)-xhat8(1))./(xhat8(3)+xhat8(4).*exp(-xhat8(5).*(indvar'-xhat8(6)))).^(1./xhat8(7))),0),1) ;
plot(indvar(1:N),proboffer_f_dat(1:N),'r',indvar(1:N),AdjustedProbabilityofReceivinganOfferWomen(1:N),'b')

[xhat8,fval8,exitflag8]         = fminsearch(@(x) adcurve6(x,indvar(1:15),(ProportionofJobSeekersActivelySearching1819cohortsMen(1:15))),[.1 .5  .5 .5 .5 .5 .6],options) ;
[probsearch_m_dat]              = min(max((xhat8(1)+(xhat8(2)-xhat8(1))./(xhat8(3)+xhat8(4).*exp(-xhat8(5).*(indvar'-xhat8(6)))).^(1./xhat8(7))),0),1) ;

[xhat8,fval8,exitflag8]         = fminsearch(@(x) adcurve6(x,indvar(1:15),(ProportionofJobSeekersActivelySearching1819cohortsWomen(1:15))),[.1 .5  .5 .5 .5 .5 .6],options) ;
[probsearch_f_dat]              = min(max((xhat8(1)+(xhat8(2)-xhat8(1))./(xhat8(3)+xhat8(4).*exp(-xhat8(5).*(indvar'-xhat8(6)))).^(1./xhat8(7))),0),1) ;
plot(indvar(1:15),probsearch_m_dat(1:15),'r',indvar(1:15),(ProportionofJobSeekersActivelySearching1819cohortsMen(1:15)),'b',...
    indvar(1:15),probsearch_f_dat(1:15),'r',indvar(1:15),(ProportionofJobSeekersActivelySearching1819cohortsWomen(1:15)),'b')




beliefs_m_dat           = [MeanBaselineExpectation(1) MeanBaselineExpectation(1)*(1+RevisioninMeanExpectMidBaseBase(1))]./scale ;
beliefs_f_dat           = [MeanBaselineExpectation(2) MeanBaselineExpectation(2)*(1+RevisioninMeanExpectMidBaseBase(2))]./scale ;


for t = 1:capT2
    jfr_m_dat(t,1) = (cumsum_m_dat(t+1)-cumsum_m_dat(t))./(1-cumsum_m_dat(t)) ;
    jfr_f_dat(t,1) = (cumsum_f_dat(t+1)-cumsum_f_dat(t))./(1-cumsum_f_dat(t)) ;
end

P           = sobolset(8) ; 
x0          = net(P,100000)' ;
miniota     = 1.1 ;
maxiota     = 2.5 ;
minmu       = .01;
maxmu       = .5 ;
mingam      = .1 ;
maxgam      = .5 ;
minz        = .01 ;
maxz        = .99 ;
maxb        = .4 ;
minb        = .1 ;
minc        = .01 ;
maxc        = 10 ;
maxl        = 1 ;
minl        = .1 ;
parpool(50)
parfor z = 1:size(x0,2)
    mu0_m   = mustar_m+log(1+x0(1,z)*(maxmu-minmu)+minmu) ;
    gam_m   = x0(2,z)*(maxgam-mingam)+mingam ;
    zeta_m  = 0 ; %x0(3,z)*(maxz-minz)+minz ;
    iota_m  = x0(3,z)*(maxiota-miniota)+miniota ;
    b_m     = (x0(4,z)*(maxb-minb)+minb)*lognstat(mustar_m,sigstar_m) ;
    scost_m = x0(5,z)*(maxc-minc)+minc ;
    
    [wstar_m,cumcomp_m,cumulsum_m,meanacceptedw_m,expectedw_m,proboffer_m,probs_m,rej_atleast1_m,rejrate_m,fdist_m,fval_m] = modelsolve_endosearch(iota_m,b_m,zeta_m,mu0_m,gam_m,mustar_m,sigstar_m,scost_m,lam_m,bet,capT,capT2,Nw,Nwlong,cumcomp_m_dat,beliefs_m_dat,proboffer_m_dat,cumsum_m_dat) ;
    
    mu0_f  = mustar_f+log(1+x0(6,z)*(maxmu-minmu)+minmu) ;
    gam_f  = x0(7,z)*(maxgam-mingam)+mingam ;
    zeta_f = 0 ; %x0(3,z)*(maxz-minz)+minz ;
    iota_f = iota_m*(1+x0(8,z)*.5) ;
    b_f    = (x0(4,z)*(maxb-minb)+minb)*lognstat(mustar_f,sigstar_f) ;
    scost_f = x0(5,z)*(maxc-minc)+minc ;
    [wstar_f,cumcomp_f,cumulsum_f,meanacceptedw_f,expectedw_f,proboffer_f,probs_f,rej_atleast1_f,rejrate_f,fdist_f,fval_f] = modelsolve_endosearch(iota_f,b_f,zeta_f,mu0_f,gam_f,mustar_f,sigstar_f,scost_f,lam_f,bet,capT,capT2,Nw,Nwlong,cumcomp_f_dat,beliefs_f_dat,proboffer_f_dat,cumsum_f_dat) ;
    f(z) = [fval_m fval_f]*[fval_m fval_f]' ;
    f(z)
end
delete(gcp('nocreate'))
eval(['save(''modeldata_final_100K_b' int2str(bval) '.mat'')']) 
   

options     = optimset('Display','iter','TolX',1e-3,'TolFun',1e-4,'MaxFunEvals',1e3,'MaxIter',1e3) ;



[fval_new,order]    = sort(f) ;
iter                = 1 ;
obj(1,1)            = fval_new(1) ;
fbest(1,1)          = fval_new(1) ;                    
while iter<10
    z_star           = order(iter) ;
    mu0_m            = mustar_m+log(1+x0(1,z_star)*(maxmu-minmu)+minmu) ;
    gam_m            = x0(2,z_star)*(maxgam-mingam)+mingam ;
    zeta_m           = 0 ; %x0(3,z_star)*(maxz-minz)+minz ;
    iota_m           = x0(3,z_star)*(maxiota-miniota)+miniota ;
    b_m              = (x0(4,z_star)*(maxb-minb)+minb)*lognstat(mustar_m,sigstar_m) ;
    scost_m          = x0(5,z_star)*(maxc-minc)+minc ;

    mu0_f            = mustar_f+log(1+x0(6,z_star)*(maxmu-minmu)+minmu) ;
    gam_f            = x0(7,z_star)*(maxgam-mingam)+mingam ;
    zeta_f           = 0 ; %x0(3,z_star)*(maxz-minz)+minz ;
    iota_f           = iota_m*(1+x0(8,z_star)*.5) ;
    b_f              = (x0(4,z_star)*(maxb-minb)+minb)*lognstat(mustar_f,sigstar_f) ;
    scost_f          = x0(5,z_star)*(maxc-minc)+minc ;

    if iter==1
        vals(iter,:)     = [ mu0_m gam_m  iota_m b_m  scost_m mu0_f gam_f iota_f ] ;
        startval(iter,:) = [ log(vals(iter,1)-mustar_m) log(vals(iter,2))  log(vals(iter,3)-1.01) ...
                                sqrt(asin(vals(iter,4)./lognstat(mustar_m,sigstar_m)))  log(vals(iter,5)) log(vals(iter,6)-mustar_f) log(vals(iter,7)-vals(iter,2)) log(vals(iter,8)-vals(iter,3)+eps) ] ;
        best             = vals(iter,:) ;
        xhat_best        = startval(iter,:) ;
    else
        vals(iter,:)     = .05.*[ mu0_m gam_m  iota_m b_m  scost_m mu0_f gam_f iota_f ]+.95.*best ;
        startval(iter,:) = [ log(vals(iter,1)-mustar_m) log(vals(iter,2))  log(vals(iter,3)-1.01) ...
        sqrt(asin(vals(iter,4)./lognstat(mustar_m,sigstar_m)))  log(vals(iter,5)) log(vals(iter,6)-mustar_f) log(vals(iter,7)-vals(iter,2)) log(vals(iter,8)-vals(iter,3)+eps) ] ;
    end
    
    [xhat,obj(1,iter),exitflag]   = fminsearch(@(x) modelobjective_endosearch(x,mustar_m,mustar_f,sigstar_m,sigstar_f,bet,capT,capT2,Nw,Nwlong,...
    cumcomp_m_dat,cumcomp_f_dat,beliefs_m_dat,beliefs_f_dat,proboffer_m_dat,proboffer_f_dat,cumsum_m_dat,cumsum_f_dat,lam_m,lam_f,zeta_m,zeta_f),startval(iter,:),options) ;

    eval(['save(''modeldata_final_100K_b' int2str(bval) '_iter' int2str(iter) '.mat'')']) 

    
        if obj(1,iter)<fbest
            fbest = obj(1,iter) ;
            xhat_best = xhat ;
            [~,mu0_m,gam_m,iota_m,b_m,b_f,zeta_m,zeta_f,scost_m,scost_f,lam_m,mu0_f,gam_f,iota_f,lam_f]  = modelobjective_endosearch(xhat_best,mustar_m,mustar_f,sigstar_m,sigstar_f,bet,capT,capT2,Nw,Nwlong,...
            cumcomp_m_dat,cumcomp_f_dat,beliefs_m_dat,beliefs_f_dat,proboffer_m_dat,proboffer_f_dat,cumsum_m_dat,cumsum_f_dat,lam_m,lam_f,zeta_m,zeta_f)
            best      = [ mu0_m gam_m  iota_m b_m  scost_m mu0_f gam_f iota_f ] ;

            [wstar_m,cumcomp_m,cumulsum_m,meanacceptedw_m,expectedw_m,proboffer_m,probs_m,rej_atleast1_m,rejrate_m,fdist_m,fval_m,wgriduse_m,ecost_m,expected_off_m]  = modelsolve_endosearch(iota_m,b_m,zeta_m,mu0_m,gam_m,mustar_m,sigstar_m,scost_m,lam_m,bet,capT,capT2,Nw,Nwlong,cumcomp_m_dat,beliefs_m_dat,proboffer_m_dat,cumsum_m_dat) ;

            [wstar_f,cumcomp_f,cumulsum_f,meanacceptedw_f,expectedw_f,proboffer_f,probs_f,rej_atleast1_f,rejrate_f,fdist_f,fval_f,wgriduse_f,ecost_f,expected_off_f]  = modelsolve_endosearch(iota_f,b_f,zeta_f,mu0_f,gam_f,mustar_f,sigstar_f,scost_f,lam_f,bet,capT,capT2,Nw,Nwlong,cumcomp_f_dat,beliefs_f_dat,proboffer_f_dat,cumsum_f_dat) ;
                
        end


    Zstar(iter,:)       = best ;
    Xhat_star(iter,:)   = xhat_best; 
    
    eval(['save(''modeldata_final_100K_b' int2str(bval) '.mat'')']) 

    if iter>1
        if sum((Zstar(iter,:)-Zstar(iter,:)).^2)==0
            iter = 10 ;
            break
        end
    end
    iter = iter+1 ;

end

[~,mu0_m,gam_m,iota_m,b_m,b_f,zeta_m,zeta_f,scost_m,scost_f,lam_m,mu0_f,gam_f,iota_f,lam_f]  = modelobjective_endosearch(Xhat_star(end,:),mustar_m,mustar_f,sigstar_m,sigstar_f,bet,capT,capT2,Nw,Nwlong,...
cumcomp_m_dat,cumcomp_f_dat,beliefs_m_dat,beliefs_f_dat,proboffer_m_dat,proboffer_f_dat,cumsum_m_dat,cumsum_f_dat,lam_m,lam_f,zeta_m,zeta_f)

[wstar_m,cumcomp_m,cumulsum_m,meanacceptedw_m,expectedw_m,proboffer_m,probs_m,rej_atleast1_m,rejrate_m,fdist_m,fval_m,wgriduse_m,ecost_m,expected_off_m,cstar_m,jfr_m]  = modelsolve_endosearch(iota_m,b_m,zeta_m,mu0_m,gam_m,mustar_m,sigstar_m,scost_m,lam_m,bet,capT,capT2*10,Nw,Nwlong,cumcomp_m_dat,beliefs_m_dat,proboffer_m_dat,cumsum_m_dat) ;

[wstar_f,cumcomp_f,cumulsum_f,meanacceptedw_f,expectedw_f,proboffer_f,probs_f,rej_atleast1_f,rejrate_f,fdist_f,fval_f,wgriduse_f,ecost_f,expected_off_f,cstar_f,jfr_f]  = modelsolve_endosearch(iota_f,b_f,zeta_f,mu0_f,gam_f,mustar_f,sigstar_f,scost_f,lam_f,bet,capT,capT2*10,Nw,Nwlong,cumcomp_f_dat,beliefs_f_dat,proboffer_f_dat,cumsum_f_dat) ;


eval(['save(''modeldata_final_100K_b' int2str(bval) '.mat'')']) 



%%

% ------------------------------------- %
%  comparative statics in iota and mu0  %
% ------------------------------------- %
clear all
run cstatics_iota.m
run cstatics_mu0.m


% ------------------------------------- %
%  estimation fit, model tables         %
% ------------------------------------- %
run model_figures.m

